subroutine jagp_local_energy_adjoints_aa_aaaa_f( &
       ni, & ! (in)    number of occupied orbitals
       na, & ! (in)    number of unoccupied orbitals
       RT, & ! (in)    unoccupied row matrix times inverse of occ pairing matrix
      RTt, & ! (in)    transpose of RT
      Jsi, & ! (in)    single index occ jastrow intermediate
      Jsa, & ! (in)    single index unocc jastrow intermediate
     Jdia, & ! (in)    double index occ-unocc jastrow intermediate
     Jdai, & ! (in)    double index unocc-occ jastrow intermediate
     Jdij, & ! (in)    double index occ-occ jastrow intermediate
     Jdab, & ! (in)    double index unocc-unocc jastrow intermediate
      Tia, & ! (in)    occ-unocc slice of the one-electron integrals
    Viajb, & ! (in)    occ-unocc-occ-unocc slice of the two-electron integrals
        E, & ! (out)   energy contribution
      RTA, & ! (out)   adjoint for RT
      Ksi, & ! (out)   adjoint for single index occ jastrow intermediate
      Ksa, & ! (out)   adjoint for single index unocc jastrow intermediate
     Kdia, & ! (out)   adjoint for double index occ-unocc jastrow intermediate
     Kdai, & ! (out)   adjoint for double index unocc-occ jastrow intermediate
     Kdij, & ! (out)   adjoint for double index occ-occ jastrow intermediate
     Kdab)   ! (out)   adjoint for double index unocc-unocc jastrow intermediate

implicit none

integer,      intent(in)    ::    ni
integer,      intent(in)    ::    na
real(kind=8), intent(in)    ::    RT(na,ni)
real(kind=8), intent(in)    ::    RTt(ni,na)
real(kind=8), intent(in)    ::    Jsi(ni)
real(kind=8), intent(in)    ::    Jsa(na)
real(kind=8), intent(in)    ::    Jdia(ni,na)
real(kind=8), intent(in)    ::    Jdai(na,ni)
real(kind=8), intent(in)    ::    Jdij(ni,ni)
real(kind=8), intent(in)    ::    Jdab(na,na)
real(kind=8), intent(in)    ::    Tia(ni,na)
real(kind=8), intent(in)    ::    Viajb(ni,na,ni,na)
real(kind=8), intent(out)   ::    E
real(kind=8), intent(out)   ::    RTA(ni,na)
real(kind=8), intent(out)   ::    Ksi(ni)
real(kind=8), intent(out)   ::    Ksa(na)
real(kind=8), intent(out)   ::    Kdia(ni,na)
real(kind=8), intent(out)   ::    Kdai(na,ni)
real(kind=8), intent(out)   ::    Kdij(ni,ni)
real(kind=8), intent(out)   ::    Kdab(na,na)

integer :: i, a, j, b

real(kind=8) :: tempQ

real(kind=8) :: yjb, yj, yaj, yab, ya, yia, yij, yib, yi, yiajb
real(kind=8) :: xjb, xj, xaj, xab, xa, xia, xij, xib, xi, xiajb

! initialize outputs to zero
E = 0.00d+00
RTA = 0.00d+00
Ksi = 0.00d+00
Ksa = 0.00d+00
Kdia = 0.00d+00
Kdai = 0.00d+00
Kdij = 0.00d+00
Kdab = 0.00d+00

! one electron integrals
! begin a-loop
do a = 1,na

  ! begin i-loop
  do i = 1,ni

    ! compute intermediates
    yia = Jsa(a) * Jdia(i,a)
    yi  = yia * Jsi(i)

    ! get jastrow-modified integral
    yiajb = Tia(i,a) * yi

    ! add contribution to the RT term's adjoint
    RTA(i,a) = RTA(i,a) + yiajb

    ! get the contribution to the local energy apart from the jastrow-modified integral
    tempQ = RTt(i,a)

    ! add local energy contribution
    E = E + tempQ * yiajb

    ! get intermediate adjoints
    xiajb = tempQ
    xi    = xiajb * Tia(i,a)
    xia   = xi    * Jsi(i)

    ! add contributions to jastrow adjoints
    Ksi(i)    = Ksi(i)    + xi  * yia
    Kdia(i,a) = Kdia(i,a) + xia * Jsa(a)
    Ksa(a)    = Ksa(a)    + xia * Jdia(i,a)

  ! end i-loop
  enddo

! end a-loop
enddo


! two electron integrals
! begin b-loop
do b = 1,na

  ! begin j-loop
  do j = 1,ni

    ! compute j-loop intermediates
    yjb = Jsa(b) * Jdia(j,b)
    yj  = yjb * Jsi(j)

    ! initialize top-level j-loop intermediate adjoint
    xj = 0.00d+00

    ! begin a-loop
    do a = 1,na

      ! compute a-loop intermediates
      yab = yj  * Jdab(a,b)
      yaj = yab * Jdai(a,j)
      ya  = yaj * Jsa(a)

      ! initialize top-level a-loop intermediate adjoint
      xa = 0.00d+00

      ! begin i-loop
      do i = 1,ni

        ! compute i-loop intermediates
        yib = ya  * Jdia(i,b)
        yij = yib * Jdij(i,j)
        yia = yij * Jdia(i,a)
        yi  = yia * Jsi(i)

        ! get jastrow-modified tei element
        yiajb = Viajb(i,a,j,b) * yi

        ! add RT adjoint coulomb contributions
        RTA(j,b) = RTA(j,b) + yiajb * RTt(i,a)

        ! add RT adjoint exchange contributions
        RTA(i,b) = RTA(i,b) - yiajb * RT(a,j)

        ! get the portion of the energy contribution other than the jastrow-modified integral
        tempQ = 0.50d+00 * ( RTt(i,a) * RTt(j,b) - RT(a,j) * RTt(i,b) )

        ! add energy contribution
        E = E + yiajb * tempQ

        ! get i-loop intermediate adjoints
        xiajb = tempQ
        xi    = xiajb * Viajb(i,a,j,b)
        xia   = xi    * Jsi(i)
        xij   = xia   * Jdia(i,a)
        xib   = xij   * Jdij(i,j)

        ! add contribution to top-level a-loop intermediate adjoint
        xa = xa + xib * Jdia(i,b)

        ! add i-loop contributions to jastrow adjoints
        Ksi(i)    = Ksi(i)    + xi  * yia
        Kdia(i,a) = Kdia(i,a) + xia * yij
        Kdij(i,j) = Kdij(i,j) + xij * yib
        Kdia(i,b) = Kdia(i,b) + xib * ya

      ! end i-loop
      enddo

      ! get remaining a-loop intermediate adjoints
      xaj = xa  * Jsa(a)
      xab = xaj * Jdai(a,j)

      ! add contributions to top-level j-loop intermediate adjoint
      xj = xj + xab * Jdab(a,b)

      ! add a-loop contributions to jastrow adjoints
      Ksa(a)    = Ksa(a)    + xa  * yaj
      Kdai(a,j) = Kdai(a,j) + xaj * yab
      Kdab(a,b) = Kdab(a,b) + xab * yj

    ! end a-loop
    enddo

    ! get remaining j-loop intermediate adjoint
    xjb = xj * Jsi(j)

    ! add j-loop and b-loop contributions to jastrow adjoints
    Ksi(j)    = Ksi(j)    + xj  * yjb
    Kdia(j,b) = Kdia(j,b) + xjb * Jsa(b)
    Ksa(b)    = Ksa(b)    + xjb * Jdia(j,b)

  ! end j-loop
  enddo

! end b-loop
enddo

end subroutine jagp_local_energy_adjoints_aa_aaaa_f




subroutine jagp_local_energy_adjoints_aabb_f( &
       ni, & ! (in)    number of occupied orbitals
       na, & ! (in)    number of unoccupied orbitals
      RTt, & ! (in)    transpose of contraction of unoccupied rows with inverse matrix
       TC, & ! (in)    contraction of unoccupied columns with inverse matrix
      RTC, & ! (in)    contraction of RT with unoccupied columns
     IPMt, & ! (in)    transpose of the inverse pairing matrix
      UPM, & ! (in)    the unoccupied pairing matrix
      Jsi, & ! (in)    single index jastrow intermediate for i
      Jsa, & ! (in)    single index jastrow intermediate for a
      Jsj, & ! (in)    single index jastrow intermediate for j
      Jsb, & ! (in)    single index jastrow intermediate for b
     Jdia, & ! (in)    double index jastrow intermediate for i and a
     Jdij, & ! (in)    double index jastrow intermediate for i and j
     Jdib, & ! (in)    double index jastrow intermediate for i and b
     Jdaj, & ! (in)    double index jastrow intermediate for a and j
     Jdab, & ! (in)    double index jastrow intermediate for a and b
     Jdjb, & ! (in)    double index jastrow intermediate for j and b
    Viajb, & ! (in)    two electron integrals in 1122 order
        E, & ! (out)   energy contribution
      VRT, & ! (out)   contraction of integrals with the RT matrix
      VTC, & ! (out)   contraction of integrals with the TC matrix
       VT, & ! (out)   contraction of integrals with inverse pairing matrix
    VARTC, & ! (out)   contraction of integrals with ( UPM - RTC )
      Ksi, & ! (out)   adjoint for single index jastrow intermediate for i
      Ksa, & ! (out)   adjoint for single index jastrow intermediate for a
      Ksj, & ! (out)   adjoint for single index jastrow intermediate for j
      Ksb, & ! (out)   adjoint for single index jastrow intermediate for b
     Kdia, & ! (out)   adjoint for double index jastrow intermediate for i and a
     Kdij, & ! (out)   adjoint for double index jastrow intermediate for i and j
     Kdib, & ! (out)   adjoint for double index jastrow intermediate for i and b
     Kdaj, & ! (out)   adjoint for double index jastrow intermediate for a and j
     Kdab, & ! (out)   adjoint for double index jastrow intermediate for a and b
     Kdjb)   ! (out)   adjoint for double index jastrow intermediate for j and b

implicit none

integer,      intent(in)    ::   ni
integer,      intent(in)    ::   na
real(kind=8), intent(in)    ::   RTt(ni,na)
real(kind=8), intent(in)    ::   TC(ni,na)
real(kind=8), intent(in)    ::   RTC(na,na)
real(kind=8), intent(in)    ::   IPMt(ni,ni)
real(kind=8), intent(in)    ::   UPM(na,na)
real(kind=8), intent(in)    ::   Jsi(ni)
real(kind=8), intent(in)    ::   Jsa(na)
real(kind=8), intent(in)    ::   Jsj(ni)
real(kind=8), intent(in)    ::   Jsb(na)
real(kind=8), intent(in)    ::   Jdia(ni,na)
real(kind=8), intent(in)    ::   Jdij(ni,ni)
real(kind=8), intent(in)    ::   Jdib(ni,na)
real(kind=8), intent(in)    ::   Jdaj(na,ni)
real(kind=8), intent(in)    ::   Jdab(na,na)
real(kind=8), intent(in)    ::   Jdjb(ni,na)
real(kind=8), intent(in)    ::   Viajb(ni,na,ni,na)
real(kind=8), intent(out)   ::   E
real(kind=8), intent(out)   ::   VRT(ni,na)
real(kind=8), intent(out)   ::   VTC(ni,na)
real(kind=8), intent(out)   ::   VT(na,na)
real(kind=8), intent(out)   ::   VARTC(ni,ni)
real(kind=8), intent(out)   ::   Ksi(ni)
real(kind=8), intent(out)   ::   Ksa(na)
real(kind=8), intent(out)   ::   Ksj(ni)
real(kind=8), intent(out)   ::   Ksb(na)
real(kind=8), intent(out)   ::   Kdia(ni,na)
real(kind=8), intent(out)   ::   Kdij(ni,ni)
real(kind=8), intent(out)   ::   Kdib(ni,na)
real(kind=8), intent(out)   ::   Kdaj(na,ni)
real(kind=8), intent(out)   ::   Kdab(na,na)
real(kind=8), intent(out)   ::   Kdjb(ni,na)

integer :: i, a, j, b

real(kind=8) :: tempQ, tempARTC

real(kind=8) :: yjb, yj, yaj, yab, ya, yia, yij, yib, yi, yiajb
real(kind=8) :: xjb, xj, xaj, xab, xa, xia, xij, xib, xi, xiajb

E = 0.00d+00
VRT = 0.00d+00
VTC = 0.00d+00
VT = 0.00d+00
VARTC = 0.00d+00
Ksi = 0.00d+00
Ksa = 0.00d+00
Ksj = 0.00d+00
Ksb = 0.00d+00
Kdia = 0.00d+00
Kdij = 0.00d+00
Kdib = 0.00d+00
Kdaj = 0.00d+00
Kdab = 0.00d+00
Kdjb = 0.00d+00

! begin b-loop
do b = 1,na

  ! begin j-loop
  do j = 1,ni

    ! compute j-loop intermediates
    yjb = Jsb(b) * Jdjb(j,b)
    yj  = yjb * Jsj(j)

    ! initialize top-level j-loop intermediate adjoint
    xj = 0.00d+00

    ! begin a-loop
    do a = 1,na

      ! compute a-loop intermediates
      yab = yj  * Jdab(a,b)
      yaj = yab * Jdaj(a,j)
      ya  = yaj * Jsa(a)
      tempARTC = UPM(a,b) - RTC(a,b)

      ! initialize top-level a-loop intermediate adjoint
      xa = 0.00d+00

      ! begin i-loop
      do i = 1,ni

        ! compute i-loop intermediates
        yib = ya  * Jdib(i,b)
        yij = yib * Jdij(i,j)
        yia = yij * Jdia(i,a)
        yi  = yia * Jsi(i)

        ! get jastrow-modified integral
        yiajb = Viajb(i,a,j,b) * yi

        ! remember the rest of the local energy term
        tempQ = RTt(i,a) * TC(j,b) + IPMt(i,j) * tempARTC

        ! process terms involving integral times RT
        VRT(j,b) = VRT(j,b) + yiajb * RTt(i,a)

        ! process terms involving integral times TC
        VTC(i,a) = VTC(i,a) + yiajb * TC(j,b)

        ! process terms involving integral times T (note that T is the inverse pairing matrix)
        VT(a,b) = VT(a,b) + yiajb * IPMt(i,j)

        ! process terms involving integral times ARTC term
        VARTC(i,j) = VARTC(i,j) + yiajb * tempARTC

        ! add energy contribution
        E = E + yiajb * tempQ

        ! get i-loop intermediate adjoints
        xiajb = tempQ
        xi    = xiajb * Viajb(i,a,j,b)
        xia   = xi    * Jsi(i)
        xij   = xia   * Jdia(i,a)
        xib   = xij   * Jdij(i,j)

        ! add contribution to top-level a-loop intermediate adjoint
        xa = xa + xib * Jdib(i,b)

        ! add i-loop contributions to jastrow adjoints
        Ksi(i)    = Ksi(i)    + xi  * yia
        Kdia(i,a) = Kdia(i,a) + xia * yij
        Kdij(i,j) = Kdij(i,j) + xij * yib
        Kdib(i,b) = Kdib(i,b) + xib * ya

      ! end i-loop
      enddo

      ! get remaining a-loop intermediate adjoints
      xaj = xa  * Jsa(a)
      xab = xaj * Jdaj(a,j)

      ! add contributions to top-level j-loop intermediate adjoint
      xj = xj + xab * Jdab(a,b)

      ! add a-loop contributions to jastrow adjoints
      Ksa(a)    = Ksa(a)    + xa  * yaj
      Kdaj(a,j) = Kdaj(a,j) + xaj * yab
      Kdab(a,b) = Kdab(a,b) + xab * yj

    ! end a-loop
    enddo

    ! get remaining j-loop intermediate adjoint
    xjb = xj * Jsj(j)

    ! add j-loop and b-loop contributions to jastrow adjoints
    Ksj(j)    = Ksj(j)    + xj  * yjb
    Kdjb(j,b) = Kdjb(j,b) + xjb * Jsb(b)
    Ksb(b)    = Ksb(b)    + xjb * Jdjb(j,b)

  ! end j-loop
  enddo

! end b-loop
enddo

end subroutine jagp_local_energy_adjoints_aabb_f
